FIE5-2-IOS-4-Devoir-2

Etude de cas sur l’analyse et la mesure de performance des Flux de Patients

Introduction

Ce travail porte sur l’analyse des flux de patients sur un plateau mutualisé de consultations externes d’un hôpital. L’objectif est de réaliser un diagnostic objectif de la performance organisationnelle du service d’urologie en s’appuyant sur les méthodes et outils d’analyse de flux vus précédemment ainsi que sur le langage R à travers l’environnement de RStudio. La figure ci-dessous illustre les principaux flux ainsi que le Plan du plateau de consultations (voir fichiers PlanConsultations.pdf et ZoomURO.pdf consultables à partir http://bit.ly/PlansCHUTlse)

Afin de collecter des données sur les parcours suivis, les patients qui se sont présentés le 12/11/2015 ont été équipés d’une étiquette électronique (type RFID) qui a permis de tracer leurs parcours dans le plateau de consultation. Les données collectées ont été fusionnées avec les données des outils de gestion des dossiers administratifs et médicaux utilisés par les personnels. L’ensemble est disponible sous la forme d’un fichier log, illustré par le tableau suivant (voir annexe LogPatientUROseul_12112015.xlsx consultable à partir http://bit.ly/logPatients).

Vue du tableau de données de log patient

Les différentes de ce tableau de données sont :

• ID (Col A) : Identifiant du patient
• Timestamp start (Col B) : horodatage entrée de zone ou salle
• Timestamp end (Col C) : horodatage sortie de zone ou salle
• Activity_MACRO (Col D) : Activité suivie et indice salle (i)
• Activity_DETAILS (Col E) : Type d’activité
• Ress.Humaines (Col F) : Ressources humaines administratives ou soignantes intervenant dans l’activité
• distance parcourues (Col G) : Distance parcourue cumulée
• début/fin opX (Col H à O) : horodatage début/fin de chaque opération de prise en charge par une ressources administrative ou soignante (maxi 4 opérations par activité).

Devoir 2 - Travail demandé : Modélisation data-driven du service d’urologie

Ce travail est à rendre pour le 11/02/2026 (avant le cours) au format Rmd ou R (+pdf si besoin) avec pour titre “NOM_Prenom_FIE5-2-IOS-4-Devoir-2.*”

Ce deuxième devoir se place en suite directe du premier et vous propose un premier type de modélisation data-driven (et top down), c’est-à-dire basé sur les travaux de Whitt & Zhang (2017) A data-driven model of an emergency department.

Question 1) Loi de Little (/3)

Pour rappel, la loi de Little est une loi fondamentale qui, dans un cadre asymptotique, lie le niveau d’occupation moyen au temps d’attente moyen par le taux d’arrivée moyen selon la formule donnée ci-dessous : \[ L = \lambda W \] - \(L\) : Niveau d’occupation moyen (asymptotique) - \(\lambda\) : Taux d’arrivée moyen (asymptotique) - \(W\) : Temps d’attente moyen (asymptotique)

Cette loi s’applique quel que soit le système ou modèle considéré à partir du moment où \(L<+\infty\), \(W<+\infty\) et \(\lambda<+\infty\) (Little JDC, Graves SC. Little’s law. In: International series in operations research and management science. 2008: 81–100.)

Pour cette première question, il vous est demandé de vérifier que cette loi “s’applique” bien sur notre service d’urologie durant la période de 8h à 18h de la journée du 12/11/2015.

Réponse :

Attention : utiliser bien les **arrivées initiales** des patients et pas les arrivées intermédiaires après une transition

• Calcul le niveau d’occupation moyen :

  Rappel et conseil : \(L = \frac{1}{\int dt} \int l(t)dt\) avec \(l(t) = A(t) - D(t)\), vous pouvez calculer en décomposant la période en zone \(i\) de même niveau d’occupation et en calculant \(L[8h;18h] = \frac{1}{\sum_i t_i} \sum_i l_i t_i\) avec \(t_i\) la durée de la zone \(i\)

library(readxl)
library(dplyr)

# Calcul de L
# Load the Excel dataset Patient URO
df <- read_excel("Log_Patient_URO_12112015.xlsx")

# Rename columns for easier handling
df1 <- df %>%
  rename(
    Ress_Humaines = `Ress. Humaines`,
    Timestamp_start = `Timestamp start`,
    Timestamp_end   = `Timestamp end`,
    DISTANCE_PARCOURUE = `distance parcourue`
  ) %>%
  mutate(
    Timestamp_start = as.POSIXct(Timestamp_start,
                                 format = "%d/%m/%Y %H:%M:%S",
                                 tz = "Europe/Paris")
  )

# Get time in system based on entry and exit
temps_systeme <- df1 %>%
  group_by(ID) %>%
  summarise(
    entree = min(Timestamp_start, na.rm = TRUE),
    sortie = max(Timestamp_start,   na.rm = TRUE)
  ) %>%
  ungroup()

events_arr <- temps_systeme %>%
  dplyr::select(time = entree) %>%
  mutate(delta = 1)

events_dep <- temps_systeme %>%
  dplyr::select(time = sortie) %>%
  mutate(delta = -1)

events <- bind_rows(events_arr, events_dep) %>%
  arrange(time)

events <- events %>%
  mutate(
    L_t = cumsum(delta),
    dt = as.numeric(difftime(lead(time), time, units = "hours"))
  )

L <- sum(events$L_t * events$dt, na.rm = TRUE) /
     sum(events$dt, na.rm = TRUE)

L

[1] 9.134184

• Calcul du taux d’arrivée moyen :

library(dplyr)
library(readxl)
library(tidyverse)
# Calcul de lamba
# Load the Excel dataset Patient URO
df <- read_excel("Log_Patient_URO_12112015.xlsx")
df

# Rename columns for easier handling
df1 <- df %>%
   rename(Ress_Humaines = `Ress. Humaines`,
        Timestamp_start = `Timestamp start`,
        Timestamp_end   = `Timestamp end`,
        DISTANCE_PARCOURUE = `distance parcourue`)

arrivees <- df1 %>%
  filter(Activity_MACRO == "Entrée des Consultations") %>%
  group_by(ID) %>%
  summarise(arrivee = min(Timestamp_start)) %>%
  ungroup()

arrivees

arrivees_8_18 <- arrivees %>%
  filter(format(arrivee, "%H:%M:%S") >= "08:00:00",
         format(arrivee, "%H:%M:%S") <= "18:00:00")

arrivees_8_18

lambda <- nrow(arrivees_8_18) / 10
lambda

[1] 6.8

• Calcul du niveau d’attente moyen :

# Load the Excel dataset Patient URO
df <- read_excel("Log_Patient_URO_12112015.xlsx")
df

# Rename columns for easier handling
df1 <- df %>%
   rename(Ress_Humaines = `Ress. Humaines`,
        Timestamp_start = `Timestamp start`,
        Timestamp_end   = `Timestamp end`,
        DISTANCE_PARCOURUE = `distance parcourue`)

# Get time in system based on entry and exit
temps_systeme <- df1 %>%
  group_by(ID) %>%
  summarise(
    entree = min(Timestamp_start, na.rm = TRUE),
    sortie = max(Timestamp_start, na.rm = TRUE)
  ) %>%
  ungroup()

temps_systeme

temps_systeme <- temps_systeme %>%
  mutate(W = as.numeric(difftime(sortie, entree, units = "hours")))

W <- mean(temps_systeme$W, na.rm = TRUE)
W

[1] 1.772818

• Comparaison de \(L\) et \(\lambda W\) :

```r
library(dplyr)
library(readxl)

# OCCUPATION
events <- bind_rows(
  temps_systeme %>% transmute(time = entree, delta = 1),
  temps_systeme %>% transmute(time = sortie, delta = -1)
) %>%
  arrange(time) %>%
  mutate(
    L_t = cumsum(delta),
    dt = as.numeric(difftime(lead(time), time, units = "hours"))
  )

L <- sum(events$L_t * events$dt, na.rm = TRUE) /
     sum(events$dt, na.rm = TRUE)

L 
```

```
[1] 9.134184
```

```r
lambda*W
```

```
[1] 12.05516
```

• Rappel et conseil : Tenter d’expliquer pourquoi on observe une différence. En particulier, pour cette explication penser à pourquoi \(\lambda W\) peut être assimilée à une moyenne des niveaux de présence moyen de patients \(m\) : \(\frac{1}{M}\sum_{m \in M} \frac{1}{t}\int_0^t \mathbb{1}_m(t) dt\) avec \(\mathbb{1}_m(t)\) la fonction d’indicatrice de présence du patient \(m\)

Réponse :

On observe que L≈9 alors que λ*W≈12 La différence s’explique par le fait que :

• Effets de bord temporels : La journée ne couvre pas un cycle complet (patients arrivés avant 8h ou partis après 18h)

• Non-stationnarité : Les arrivées sont concentrées le matin, les départs étalés

• Moyenne temporelle vs moyenne par patient : L est une moyenne dans le temps, λW est une moyenne sur les patients

Question 2) Processus de Poisson d’arrivée non homogène (/3)

Un processus de poisson est un processus de comptage (dans le temps) indiquant un nombre évènements ayant occurés entre un temps \(0\) et un temps \(t\) selon une distribution de Poisson \(\mathcal(P)(\lambda * t)\) avec un taux par unité de temps \(\lambda\). Nous verrons dans le cours 3 que les temps entre chaque événement suive une loi de distribution exponentielle de paramètre \(\lambda\).

Un processus de Poisson non homogène (NHPP) est un processus de comptage où le taux d’événement \(\lambda(t)\) n’est pas constant. En considérant une modélisation de NHPP basé sur taux d’arrivées des tranches [8h;10h], ]10h;12h], …, ]16h;18h] du service d’urologie générer 30 échantillons de ce NHPP sur la période 8h, 18h et illustrés les.

Pour vous aider voici un exemple pour un Processus de Poisson homogène : (n’hésitez pas à aller plus loin aussi pour l’illustration)

lambda = 10 # par heure
samples <- tibble(run=to_vec(for(i in 1:10) rep(i,200)),
       id=rep(1:200,10),
       delta_t=rexp(2000,lambda)) %>%
  group_by(run) %>%
  mutate(t = cumsum(delta_t)) %>%
  filter(t <= 10)

# illustration 1
samples %>% 
  ggplot(aes(t,id,color=factor(run))) +
  geom_point() 

# illustration 2
samples %>%   
  mutate(t = cut(t,seq(0,30,by=1),include.lowest = TRUE)) %>%
  count(run,t) %>%
  arrange(run,t) %>%
  ggplot(aes(t,n)) +
  geom_boxplot()

Attention : pour un processus de Poisson non homogène, il faudra générer le nombre d'arrivées avec _rpoiss_ pour chaque tranche de temps avec un taux différent et ensuite les répartir uniformément dans le temps (au sein de leur intervalle) en générant des valuers uniform (_runif_).

Réponse :

# Import des données
data_raw <- read_excel("Log_Patient_URO_12112015.xlsx")

# on garde que l'arrivée initiale
data_arrivals <- data_raw %>%
  mutate(arrival_time = ymd_hms(`Timestamp start`)) %>%
  group_by(ID) %>%                
  slice_min(arrival_time, n = 1) %>%        
  ungroup() %>%
  mutate(
    hour = hour(arrival_time) + minute(arrival_time)/60
  ) %>%
  filter(hour >= 8, hour < 18)

# Découpage en tranches horaires
data_arrivals <- data_arrivals %>%
  mutate(
    tranche = cut(
      hour,
      breaks = c(8,10,12,14,16,18),
      include.lowest = TRUE,
      right = TRUE
    )
  )

# Estimation des taux λ(t)
lambda_hat <- data_arrivals %>%
  count(tranche) %>%
  mutate(lambda = n / 2) 

# Définition des intervalles
intervals <- tibble(
  start  = c(8,10,12,14,16),
  end    = c(10,12,14,16,18),
  lambda = lambda_hat$lambda
)

# Simulation du NHPP (30 runs)
set.seed(123)
n_runs <- 30

samples_nhpp <- map_dfr(1:n_runs, function(run_id) {

  map_dfr(1:nrow(intervals), function(i) {

    dt <- intervals$end[i] - intervals$start[i]
    n_events <- rpois(1, intervals$lambda[i] * dt)

    tibble(
      run = run_id,
      t = runif(n_events, intervals$start[i], intervals$end[i])
    )
  })
}) %>%
  arrange(run, t) %>%
  group_by(run) %>%
  mutate(id = row_number()) %>%
  ungroup()

# Illustration 1 : trajectoires
samples_nhpp %>%
  ggplot(aes(x = t, color = factor(run))) +
  stat_ecdf(geom = "step", alpha = 0.5) +
  scale_y_continuous(labels = scales::percent) +
  labs(
    title = "Trajectoires du NHPP – Fonction de comptage",
    x = "Temps (heures)",
    y = "Proportion d'arrivées cumulées"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

# Calcul des arrivées réelles par heure
real_arrivals <- data_arrivals %>%
  mutate(hour = cut(hour, seq(8,18,by=1), include.lowest = TRUE)) %>%
  count(hour)

# Illustration 2 AMÉLIORÉE
samples_nhpp %>%
  mutate(hour = cut(t, seq(8,18,by=1), include.lowest = TRUE)) %>%
  count(run, hour) %>%
  ggplot(aes(hour, n)) +
  geom_boxplot(fill = "lightblue", alpha = 0.6) +
  geom_point(data = real_arrivals, aes(hour, n), 
             color = "red", size = 3, shape = 18) +  
  labs(
    title = "Distribution des arrivées par heure",
    subtitle = "Boxplot: NHPP simulé (30 runs) | Losanges rouges: données réelles",
    x = "Heure",
    y = "Nombre d'arrivées"
  ) +
  theme_minimal()

Question 3) Modélisation de la durée de séjour par une distribution (/3)

A l’aide de la documentation “Ricci-distributions-en.pdf” fournie, notamment avec la fonction fitdistr() de la librairie MASS, tester la modélisation (par MLE, Maximum Likelihood Estimation) la durée de séjours de l’ensemble des patients avec différentes distributions : normale, exponentielle, gamma, weibull.

Réponse :

Entrez votre texte ici

#_Entrez votre code R ici_

library(readxl)
library(MASS)

data_service <- read_excel("Log_Patient_URO_12112015.xlsx")

data_raw <- data_service %>%
  mutate(
    datetime_begin = ymd_hms(`Timestamp start`),
    datetime_end   = ymd_hms(`Timestamp end`)
  )

data_patient <- data_raw %>%
  group_by(ID) %>%
  summarise(
    arrival_time = min(datetime_begin, na.rm = TRUE),
    departure_time = max(datetime_end, na.rm = TRUE),
    .groups = "drop"
  )

data_patient <- data_patient %>%
  mutate(
    duree_sejour = as.numeric(difftime(departure_time,
                                       arrival_time,
                                       units = "hours"))
  )

duree <- data_patient$duree_sejour
duree <- duree[!is.na(duree)]
duree <- duree[duree > 0]

summary(duree)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.4654  0.7964  1.5256  1.7322  2.4347  5.8699 

hist(duree, breaks = 30)

fit_norm <- fitdistr(duree, "normal")
fit_norm

      mean          sd    
  1.73217947   1.12042754 
 (0.13587180) (0.09607587)

fit_exp <- fitdistr(duree, "exponential")
fit_exp

     rate   
  0.5773074 
 (0.0700088)

fit_gamma <- fitdistr(duree, "gamma")
fit_gamma

     shape       rate   
  2.6808643   1.5476811 
 (0.4342089) (0.2756369)

fit_weib <- fitdistr(duree, "weibull")
fit_weib

     shape       scale  
  1.6654502   1.9525277 
 (0.1507457) (0.1506968)

AIC_fitdistr <- function(fit) {
  k <- length(fit$estimate)
  -2 * fit$loglik + 2 * k
}

AIC_norm  <- AIC_fitdistr(fit_norm)
AIC_exp   <- AIC_fitdistr(fit_exp)
AIC_gamma <- AIC_fitdistr(fit_gamma)
AIC_weib  <- AIC_fitdistr(fit_weib)

AIC_values <- data.frame(
  Distribution = c("Normale", "Exponentielle", "Gamma", "Weibull"),
  AIC = c(AIC_norm, AIC_exp, AIC_gamma, AIC_weib)
)

AIC_values[order(AIC_values$AIC), ]

hist(duree, prob = TRUE, breaks = 30,
     main = "Ajustement des lois – Durée de séjour",
     xlab = "Durée")

curve(dnorm(x,
            mean = fit_norm$estimate[1],
            sd   = fit_norm$estimate[2]),
      add = TRUE, col = "blue", lwd = 2)

curve(dexp(x,
           rate = fit_exp$estimate),
      add = TRUE, col = "red", lwd = 2)

curve(dgamma(x,
             shape = fit_gamma$estimate["shape"],
             rate  = fit_gamma$estimate["rate"]),
      add = TRUE, col = "green", lwd = 2)

curve(dweibull(x,
               shape = fit_weib$estimate["shape"],
               scale = fit_weib$estimate["scale"]),
      add = TRUE, col = "purple", lwd = 2)

legend("topright",
       legend = c("Normale", "Exponentielle", "Gamma", "Weibull"),
       col = c("blue", "red", "green", "purple"),
       lwd = 2)

Quelles est la meilleure distribution ? (Justifiez) Comparez aussi les moyennes, écart-types et coefficients de variation obtenus pour chaque distribution et par rapport au calcul direct des indicateurs sur les variables statistiques.

Réponse :

Parmi les distributions étudiées, la loi de Weibull apparaît comme la plus appropriée pour modéliser la durée de séjour. Elle respecte le support strictement positif de la variable et reproduit correctement les principaux indicateurs statistiques, en particulier la variabilité observée. La loi Gamma aurait également pu constituer un choix pertinent, car elle présente des caractéristiques proches et un bon accord avec les données empiriques ; toutefois, la Weibull offre une flexibilité légèrement supérieure, notamment dans la modélisation des durées extrêmes et dans l’interprétation du comportement du taux de sortie au cours du séjour. La loi normale, bien que numériquement proche, reste conceptuellement inadaptée en raison de son support non borné inférieur, tandis que la loi exponentielle est clairement inappropriée car elle impose une variabilité trop élevée par rapport aux observations. Ainsi, la loi de Weibull est retenue comme meilleur compromis, la loi Gamma pouvant être considérée comme une alternative valable en second choix.

#_Entrez votre code R ici_

# Indicateurs empiriques
mean_emp <- mean(duree)
sd_emp   <- sd(duree)
cv_emp   <- sd_emp / mean_emp

mean_norm <- fit_norm$estimate["mean"]
sd_norm   <- fit_norm$estimate["sd"]
cv_norm   <- sd_norm / mean_norm

rate_exp <- fit_exp$estimate["rate"]

mean_exp <- 1 / rate_exp
sd_exp   <- 1 / rate_exp
cv_exp   <- sd_exp / mean_exp

shape_g <- fit_gamma$estimate["shape"]
rate_g  <- fit_gamma$estimate["rate"]

mean_gamma <- shape_g / rate_g
sd_gamma   <- sqrt(shape_g) / rate_g
cv_gamma   <- sd_gamma / mean_gamma

shape_w <- fit_weib$estimate["shape"]
scale_w <- fit_weib$estimate["scale"]

mean_weib <- scale_w * gamma(1 + 1 / shape_w)
sd_weib   <- scale_w * sqrt(
  gamma(1 + 2 / shape_w) - gamma(1 + 1 / shape_w)^2
)
cv_weib <- sd_weib / mean_weib

comparaison_table <- data.frame(
  Distribution = c("Empirique", "Normale", "Exponentielle", "Gamma", "Weibull"),
  Moyenne = c(mean_emp, mean_norm, mean_exp, mean_gamma, mean_weib),
  Ecart_type = c(sd_emp, sd_norm, sd_exp, sd_gamma, sd_weib),
  Coefficient_variation = c(cv_emp, cv_norm, cv_exp, cv_gamma, cv_weib)
)

comparaison_table

NA

Quel que soit votre réponse, refaites ce travail pour la distribution gamma en séparant la modélisation des durées de séjour des patients prioritaires et non prioritaires et comparez les deux distributions graphiquement et à l’aide leur moments (espérance, variance, coefficient de variation) et/ou leur paramètres d’échelle et de forme.

Réponse :

L’analyse des durées de séjour des patients selon une loi Gamma montre que, bien que la durée moyenne soit très proche entre les patients prioritaires (≈1,76) et non prioritaires (≈1,72), la dispersion diffère sensiblement. Les patients prioritaires présentent une variance plus élevée (1,66 contre 0,91) et un coefficient de variation plus important (0,73 contre 0,55), indiquant une plus grande variabilité relative de leurs durées de séjour. Cette différence se reflète également dans les paramètres de la loi Gamma : le paramètre shape des patients prioritaires (≈1,86) est plus faible que celui des non prioritaires (≈3,26), ce qui traduit une distribution plus asymétrique et étalée, avec une queue plus longue à droite. À l’inverse, la distribution des non prioritaires est plus concentrée autour de la moyenne. Ainsi, même si les durées moyennes sont similaires, les patients prioritaires montrent des séjours plus hétérogènes, ce qui peut avoir des implications pour la planification hospitalière et la gestion des ressources.

#_Entrez votre code R ici_

data_patient <- data_raw %>%
  mutate(
    datetime_begin = ymd_hms(`Timestamp start`),
    datetime_end   = ymd_hms(`Timestamp end`)
  ) %>%
  group_by(ID) %>%
  summarise(
    arrival_time = min(datetime_begin, na.rm = TRUE),
    departure_time = max(datetime_end, na.rm = TRUE),
    prioritaire = ifelse(any(grepl("PRIO", Activity_DETAILS)), "Oui", "Non"),
    .groups = "drop"
  ) %>%
  mutate(
    duree_sejour = as.numeric(difftime(departure_time,
                                       arrival_time,
                                       units = "hours"))
  )

duree_prio <- data_patient %>%
  filter(prioritaire == "Oui") %>%
  pull(duree_sejour)

duree_non_prio <- data_patient %>%
  filter(prioritaire == "Non") %>%
  pull(duree_sejour)

# Nettoyage
duree_prio <- duree_prio[duree_prio > 0 & !is.na(duree_prio)]
duree_non_prio <- duree_non_prio[duree_non_prio > 0 & !is.na(duree_non_prio)]

library(MASS)

fit_gamma_prio <- fitdistr(duree_prio, "gamma")
fit_gamma_non_prio <- fitdistr(duree_non_prio, "gamma")

params_gamma <- data.frame(
  Groupe = c("Prioritaires", "Non prioritaires"),
  Shape = c(fit_gamma_prio$estimate["shape"],
            fit_gamma_non_prio$estimate["shape"]),
  Rate = c(fit_gamma_prio$estimate["rate"],
           fit_gamma_non_prio$estimate["rate"])
)

params_gamma

moments_gamma <- data.frame(
  Groupe = c("Prioritaires", "Non prioritaires"),
  Esperance = c(
    fit_gamma_prio$estimate["shape"] / fit_gamma_prio$estimate["rate"],
    fit_gamma_non_prio$estimate["shape"] / fit_gamma_non_prio$estimate["rate"]
  ),
  Variance = c(
    fit_gamma_prio$estimate["shape"] / fit_gamma_prio$estimate["rate"]^2,
    fit_gamma_non_prio$estimate["shape"] / fit_gamma_non_prio$estimate["rate"]^2
  ),
  Coefficient_variation = c(
    1 / sqrt(fit_gamma_prio$estimate["shape"]),
    1 / sqrt(fit_gamma_non_prio$estimate["shape"])
  )
)

moments_gamma

hist(duree_prio, prob = TRUE, breaks = 30,
     col = rgb(1,0,0,0.35),
     xlim = range(c(duree_prio, duree_non_prio)),
     main = "Comparaison des durées de séjour – Loi Gamma",
     xlab = "Durée (heures)")

curve(dgamma(x,
             shape = fit_gamma_prio$estimate["shape"],
             rate  = fit_gamma_prio$estimate["rate"]),
      col = "red", lwd = 2, add = TRUE)

hist(duree_non_prio, prob = TRUE, breaks = 30,
     col = rgb(0,0,1,0.35),
     add = TRUE)

curve(dgamma(x,
             shape = fit_gamma_non_prio$estimate["shape"],
             rate  = fit_gamma_non_prio$estimate["rate"]),
      col = "blue", lwd = 2, add = TRUE)

legend("topright",
       legend = c("Prioritaires", "Non prioritaires"),
       col = c("red", "blue"),
       lwd = 2)

# Q-Q plot pour prioritaires
par(mfrow = c(1, 2))

# Prioritaires
qqplot(qgamma(ppoints(length(duree_prio)),
              shape = fit_gamma_prio$estimate["shape"],
              rate = fit_gamma_prio$estimate["rate"]),
       duree_prio,
       main = "Q-Q Plot - Prioritaires",
       xlab = "Quantiles théoriques (Gamma)",
       ylab = "Quantiles observés")
abline(0, 1, col = "red")

# Non-prioritaires
qqplot(qgamma(ppoints(length(duree_non_prio)),
              shape = fit_gamma_non_prio$estimate["shape"],
              rate = fit_gamma_non_prio$estimate["rate"]),
       duree_non_prio,
       main = "Q-Q Plot - Non prioritaires",
       xlab = "Quantiles théoriques (Gamma)",
       ylab = "Quantiles observés")
abline(0, 1, col = "blue")

par(mfrow = c(1, 1))

Question 4) Modélisation de la durée de séjour par régression linéaire (/3)

A l’aide de la fonction lm, construisez et analysez un modèle de régression linéaire estimant la durée de séjour qui prennent en variable d’entrée le bloc de 2h [8h;10h], …, ]16h;18h].

Voici un petit exemple d’utilisation :

Y = c(1, 1.5, 2, 3)
X = c("A","A","B","B")

model <- lm(Y ~ X)
summary(model) 

Call:
lm(formula = Y ~ X)

Residuals:
    1     2     3     4 
-0.25  0.25 -0.50  0.50 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   1.2500     0.3953   3.162   0.0871 .
XB            1.2500     0.5590   2.236   0.1548  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.559 on 2 degrees of freedom
Multiple R-squared:  0.7143,    Adjusted R-squared:  0.5714 
F-statistic:     5 on 1 and 2 DF,  p-value: 0.1548

predict(model) # predict(model,newdata=[new_dataframe]) if new data

   1    2    3    4 
1.25 1.25 2.50 2.50 

Réponse :

Entrez votre texte ici

library(dplyr)
library(lubridate)
library(ggplot2)
library(readxl)
library(broom)  # ✅ IMPORTANT : Charger broom au début !


# Chargement des données
df <- read_excel("Log_Patient_URO_12112015.xlsx")

# Renommer les colonnes
df1 <- df %>%
  rename(
    Ress_Humaines = `Ress. Humaines`,
    Timestamp_start = `Timestamp start`,
    Timestamp_end   = `Timestamp end`,
    DISTANCE_PARCOURUE = `distance parcourue`
  ) %>%
  mutate(
    Timestamp_start = as.POSIXct(Timestamp_start, format = "%d/%m/%Y %H:%M:%S", tz = "Europe/Paris"),
    Timestamp_end   = as.POSIXct(Timestamp_end, format = "%d/%m/%Y %H:%M:%S", tz = "Europe/Paris")
  )

# Calculer la durée de séjour et l'heure d'arrivée pour chaque patient
temps_systeme <- df1 %>%
  group_by(ID) %>%
  summarise(
    entree = min(Timestamp_start, na.rm = TRUE),
    sortie = max(Timestamp_end, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    W = as.numeric(difftime(sortie, entree, units = "hours")),
    heure_entree = hour(entree) + minute(entree)/60
  )

# Créer la variable bloc_2h selon l'heure d'arrivée
temps_systeme <- temps_systeme %>%
  mutate(
    bloc_2h = case_when(
      heure_entree >= 8 & heure_entree < 10 ~ "[08h;10h]",
      heure_entree >= 10 & heure_entree < 12 ~ "[10h;12h]",
      heure_entree >= 12 & heure_entree < 14 ~ "[12h;14h]",
      heure_entree >= 14 & heure_entree < 16 ~ "[14h;16h]",
      heure_entree >= 16 & heure_entree < 18 ~ "[16h;18h]",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(bloc_2h), W > 0, !is.na(W))

temps_systeme <- temps_systeme %>%
  mutate(
    bloc_2h = factor(bloc_2h, 
                     levels = c("[08h;10h]", "[10h;12h]", "[12h;14h]", 
                                "[14h;16h]", "[16h;18h]"))
  )

# Vérification
head(temps_systeme)

summary(temps_systeme$W)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.4654  0.7963  1.5256  1.7189  2.4347  5.8699 

# ==============================================================================
# MODÈLE 1
# ==============================================================================

model1 <- lm(W ~ bloc_2h, data = temps_systeme)
summary(model1)

Call:
lm(formula = W ~ bloc_2h, data = temps_systeme)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5650 -0.8471 -0.1567  0.7813  3.7957 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)        0.9368     0.4557   2.056   0.0443 *
bloc_2h[10h;12h]   0.6305     0.5196   1.213   0.2299  
bloc_2h[12h;14h]   0.7602     0.5765   1.319   0.1923  
bloc_2h[14h;16h]   1.1374     0.5582   2.038   0.0461 *
bloc_2h[16h;18h]   1.0119     0.5344   1.894   0.0632 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.116 on 59 degrees of freedom
Multiple R-squared:  0.08116,   Adjusted R-squared:  0.01887 
F-statistic: 1.303 on 4 and 59 DF,  p-value: 0.2795

# Prédictions
temps_systeme$pred_W1 <- predict(model1)

# Visualisation
ggplot(temps_systeme, aes(x = bloc_2h, y = W)) +
  geom_jitter(width = 0.2, height = 0, alpha = 0.5, color = "lightblue", size = 2) +
  geom_point(aes(y = pred_W1), color = "red", size = 4, shape = 17) +
  stat_summary(fun = mean, geom = "point", color = "darkgreen", size = 4, shape = 18) +
  labs(
    title = "Modèle 1 : Durée de séjour selon le bloc horaire",
    subtitle = "Triangles rouges = prédictions | Losanges verts = moyennes observées",
    x = "Bloc horaire d'arrivée",
    y = "Durée de séjour W (heures)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    axis.text.x = element_text(angle = 0)
  )

# ✅ CORRECTION : Diagnostics UN PAR UN au lieu de 4 en même temps
plot(model1, which = 1, main = "Modèle 1 - Résidus vs Valeurs ajustées")

plot(model1, which = 2, main = "Modèle 1 - Q-Q plot")

plot(model1, which = 3, main = "Modèle 1 - Scale-Location")

plot(model1, which = 5, main = "Modèle 1 - Résidus vs Levier")

# ==============================================================================
# MODÈLE 2
# ==============================================================================

# Identifier les patients prioritaires
temps_systeme2 <- temps_systeme %>%
  left_join(
    df1 %>%
      group_by(ID) %>%
      summarise(prioritaire = ifelse(any(grepl("PRIO", Activity_DETAILS)), "Oui", "Non"),
                .groups = "drop"),
    by = "ID"
  )

temps_systeme2 <- temps_systeme2 %>%
  mutate(prioritaire = factor(prioritaire, levels = c("Non", "Oui")))

# Vérification
table(temps_systeme2$prioritaire)

Non Oui 
 45  19 

# Modèle 2
model2 <- lm(W ~ bloc_2h + prioritaire, data = temps_systeme2)
summary(model2)

Call:
lm(formula = W ~ bloc_2h + prioritaire, data = temps_systeme2)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5567 -0.8390 -0.1515  0.7885  3.7790 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)       0.93268    0.46253   2.016   0.0484 *
bloc_2h[10h;12h]  0.62588    0.52714   1.187   0.2399  
bloc_2h[12h;14h]  0.75690    0.58286   1.299   0.1992  
bloc_2h[14h;16h]  1.13320    0.56530   2.005   0.0497 *
bloc_2h[16h;18h]  1.00986    0.53958   1.872   0.0663 .
prioritaireOui    0.02496    0.31047   0.080   0.9362  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.126 on 58 degrees of freedom
Multiple R-squared:  0.08126,   Adjusted R-squared:  0.002062 
F-statistic: 1.026 on 5 and 58 DF,  p-value: 0.4109

# Prédictions
temps_systeme2$pred_W2 <- predict(model2)

# Visualisation
ggplot(temps_systeme2, aes(x = bloc_2h, y = W, color = prioritaire)) +
  geom_jitter(width = 0.2, height = 0, alpha = 0.4, size = 2) +
  geom_point(aes(y = pred_W2), shape = 17, size = 4) +
  geom_line(aes(y = pred_W2, group = prioritaire), linewidth = 1.2) +
  labs(
    title = "Modèle 2 : W ~ bloc_2h + prioritaire",
    subtitle = "Triangles = prédictions | Lignes parallèles montrent l'effet additif",
    x = "Bloc horaire d'arrivée",
    y = "Durée de séjour W (heures)",
    color = "Prioritaire"
  ) +
  scale_color_manual(values = c("Non" = "steelblue", "Oui" = "tomato")) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    legend.position = "top"
  )

# Diagnostics UN PAR UN
plot(model2, which = 1, main = "Modèle 2 - Résidus vs Valeurs ajustées")

plot(model2, which = 2, main = "Modèle 2 - Q-Q plot")

plot(model2, which = 3, main = "Modèle 2 - Scale-Location")

plot(model2, which = 5, main = "Modèle 2 - Résidus vs Levier")

# ==============================================================================
# COMPARAISON DES MODÈLES
# ==============================================================================

models_comparison <- tibble(
  Modèle = c("Modèle 1: W ~ bloc_2h", 
             "Modèle 2: W ~ bloc_2h + prioritaire"),
  R_squared = c(summary(model1)$r.squared, 
                summary(model2)$r.squared),
  R_squared_adj = c(summary(model1)$adj.r.squared, 
                    summary(model2)$adj.r.squared),
  AIC = c(AIC(model1), AIC(model2)),
  BIC = c(BIC(model1), BIC(model2))
)

print(models_comparison)

# A tibble: 2 × 5
  Modèle                              R_squared R_squared_adj   AIC   BIC
  <chr>                                   <dbl>         <dbl> <dbl> <dbl>
1 Modèle 1: W ~ bloc_2h                  0.0812       0.0189   203.  215.
2 Modèle 2: W ~ bloc_2h + prioritaire    0.0813       0.00206  204.  220.

# Test de comparaison (ANOVA)
anova(model1, model2)

Analysis of Variance Table

Model 1: W ~ bloc_2h
Model 2: W ~ bloc_2h + prioritaire
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     59 73.526                           
2     58 73.518  1 0.0081909 0.0065 0.9362

# ==============================================================================
# MODÈLE 3 (avec interaction)
# ==============================================================================

model3 <- lm(W ~ bloc_2h * prioritaire, data = temps_systeme2)
summary(model3)

Call:
lm(formula = W ~ bloc_2h * prioritaire, data = temps_systeme2)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.8812 -0.8660 -0.1392  0.7573  3.1781 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)  
(Intercept)                      0.99721    0.51168   1.949   0.0565 .
bloc_2h[10h;12h]                 0.67149    0.60209   1.115   0.2697  
bloc_2h[12h;14h]                 0.73947    0.66994   1.104   0.2746  
bloc_2h[14h;16h]                 0.76822    0.65226   1.178   0.2441  
bloc_2h[16h;18h]                 0.97194    0.60902   1.596   0.1163  
prioritaireOui                  -0.36222    1.25335  -0.289   0.7737  
bloc_2h[10h;12h]:prioritaireOui  0.07247    1.36330   0.053   0.9578  
bloc_2h[12h;14h]:prioritaireOui  0.23016    1.48130   0.155   0.8771  
bloc_2h[14h;16h]:prioritaireOui  1.28851    1.43589   0.897   0.3735  
bloc_2h[16h;18h]:prioritaireOui  0.28076    1.41677   0.198   0.8437  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.144 on 54 degrees of freedom
Multiple R-squared:  0.1166,    Adjusted R-squared:  -0.03063 
F-statistic: 0.7919 on 9 and 54 DF,  p-value: 0.6249

# Comparaison des 3 modèles
anova(model1, model2, model3)

Analysis of Variance Table

Model 1: W ~ bloc_2h
Model 2: W ~ bloc_2h + prioritaire
Model 3: W ~ bloc_2h * prioritaire
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     59 73.526                           
2     58 73.518  1   0.00819 0.0063 0.9372
3     54 70.690  4   2.82771 0.5400 0.7070

# Visualisation du modèle 3
temps_systeme2$pred_W3 <- predict(model3)

ggplot(temps_systeme2, aes(x = bloc_2h, y = W, color = prioritaire)) +
  geom_jitter(width = 0.2, height = 0, alpha = 0.4, size = 2) +
  geom_point(aes(y = pred_W3), shape = 17, size = 4) +
  geom_line(aes(y = pred_W3, group = prioritaire), linewidth = 1.2) +  
  labs(
    title = "Modèle 3 : W ~ bloc_2h * prioritaire (avec interaction)",
    subtitle = "Les lignes NON parallèles indiqueraient une interaction",
    x = "Bloc horaire d'arrivée",
    y = "Durée de séjour W (heures)",
    color = "Prioritaire"
  ) +
  scale_color_manual(values = c("Non" = "steelblue", "Oui" = "tomato")) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    legend.position = "top"
  )

# ==============================================================================
# ✅ CORRECTION : Tableau des coefficients (besoin de library(broom))
# ==============================================================================

coef_table <- data.frame(
  Modèle = c(rep("Modèle 1", length(coef(model1))),
             rep("Modèle 2", length(coef(model2)))),
  Terme = c(names(coef(model1)), names(coef(model2))),
  Estimate = round(c(coef(model1), coef(model2)), 3),
  Std_Error = round(c(summary(model1)$coefficients[, "Std. Error"],
                      summary(model2)$coefficients[, "Std. Error"]), 3),
  P_value = round(c(summary(model1)$coefficients[, "Pr(>|t|)"],
                    summary(model2)$coefficients[, "Pr(>|t|)"]), 4)
)

coef_table$Significatif <- ifelse(coef_table$P_value < 0.05, "***", "")

print(coef_table)

     Modèle            Terme Estimate Std_Error P_value Significatif
1  Modèle 1      (Intercept)    0.937     0.456  0.0443          ***
2  Modèle 1 bloc_2h[10h;12h]    0.630     0.520  0.2299             
3  Modèle 1 bloc_2h[12h;14h]    0.760     0.576  0.1923             
4  Modèle 1 bloc_2h[14h;16h]    1.137     0.558  0.0461          ***
5  Modèle 1 bloc_2h[16h;18h]    1.012     0.534  0.0632             
6  Modèle 2      (Intercept)    0.933     0.463  0.0484          ***
7  Modèle 2 bloc_2h[10h;12h]    0.626     0.527  0.2399             
8  Modèle 2 bloc_2h[12h;14h]    0.757     0.583  0.1992             
9  Modèle 2 bloc_2h[14h;16h]    1.133     0.565  0.0497          ***
10 Modèle 2 bloc_2h[16h;18h]    1.010     0.540  0.0663             
11 Modèle 2   prioritaireOui    0.025     0.310  0.9362             

# ==============================================================================
# INTERPRÉTATION
# ==============================================================================

cat("\n=== INTERPRÉTATION DES RÉSULTATS ===\n\n")

=== INTERPRÉTATION DES RÉSULTATS ===

# Modèle 1
cat("MODÈLE 1 (W ~ bloc_2h):\n")

MODÈLE 1 (W ~ bloc_2h):

cat(sprintf("- R² = %.3f (seulement %.1f%% de la variance expliquée)\n", 
            summary(model1)$r.squared, summary(model1)$r.squared * 100))

- R² = 0.081 (seulement 8.1% de la variance expliquée)

cat("- Aucun bloc horaire n'est significatif (p > 0.05)\n")

- Aucun bloc horaire n'est significatif (p > 0.05)

cat("- Conclusion: L'heure d'arrivée n'influence PAS significativement la durée de séjour\n\n")

- Conclusion: L'heure d'arrivée n'influence PAS significativement la durée de séjour

# Modèle 2
cat("MODÈLE 2 (W ~ bloc_2h + prioritaire):\n")

MODÈLE 2 (W ~ bloc_2h + prioritaire):

cat(sprintf("- R² = %.3f (%.1f%% de la variance expliquée)\n", 
            summary(model2)$r.squared, summary(model2)$r.squared * 100))

- R² = 0.081 (8.1% de la variance expliquée)

coef_prio <- coef(model2)["prioritaireOui"]
cat(sprintf("- Effet prioritaire: +%.2f heures (soit ~%.0f minutes)\n", 
            coef_prio, coef_prio * 60))

- Effet prioritaire: +0.02 heures (soit ~1 minutes)

# Test de significativité
p_prio <- summary(model2)$coefficients["prioritaireOui", "Pr(>|t|)"]
if (p_prio < 0.05) {
  cat(sprintf("- Cet effet est SIGNIFICATIF (p = %.4f)\n", p_prio))
} else {
  cat(sprintf("- Cet effet n'est pas significatif (p = %.4f)\n", p_prio))
}

- Cet effet n'est pas significatif (p = 0.9362)

cat("\n- Conclusion: Les patients prioritaires restent significativement plus longtemps,\n")

- Conclusion: Les patients prioritaires restent significativement plus longtemps,

cat("  probablement en raison de la complexité de leurs cas (examens complémentaires).\n")

  probablement en raison de la complexité de leurs cas (examens complémentaires).

library(dplyr)
library(lubridate)

# Crée la variable bloc 2h selon l'heure d'arrivée
temps_systeme <- temps_systeme %>%
  mutate(
    heure_entree = hour(entree),
    bloc_2h = case_when(
      heure_entree >= 8 & heure_entree < 10 ~ "[08h;10h]",
      heure_entree >= 10 & heure_entree < 12 ~ "[10h;12h]",
      heure_entree >= 12 & heure_entree < 14 ~ "[12h;14h]",
      heure_entree >= 14 & heure_entree < 16 ~ "[14h;16h]",
      heure_entree >= 16 & heure_entree < 18 ~ "[16h;18h]",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(bloc_2h))  

temps_systeme

model1 <- lm(W ~ bloc_2h, data = temps_systeme)
summary(model1)

Call:
lm(formula = W ~ bloc_2h, data = temps_systeme)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5650 -0.8471 -0.1567  0.7813  3.7957 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)        0.9368     0.4557   2.056   0.0443 *
bloc_2h[10h;12h]   0.6305     0.5196   1.213   0.2299  
bloc_2h[12h;14h]   0.7602     0.5765   1.319   0.1923  
bloc_2h[14h;16h]   1.1374     0.5582   2.038   0.0461 *
bloc_2h[16h;18h]   1.0119     0.5344   1.894   0.0632 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.116 on 59 degrees of freedom
Multiple R-squared:  0.08116,   Adjusted R-squared:  0.01887 
F-statistic: 1.303 on 4 and 59 DF,  p-value: 0.2795

# prédictions sur le même dataset
temps_systeme$pred_W1 <- predict(model1)

library(ggplot2)

ggplot(temps_systeme, aes(x = bloc_2h, y = W)) +
  geom_jitter(width = 0.2, height = 0, alpha = 0.5, color = "blue") +
  geom_point(aes(y = pred_W1), color = "red", size = 3) +
  stat_summary(fun = mean, geom = "point", color = "darkgreen", size = 3) +
  labs(
    title = "Durée de séjour selon le bloc 2h",
    x = "Bloc horaire d'arrivée",
    y = "Durée de séjour W (heures)"
  ) +
  theme_minimal()

Refaites un deuxième modèle linéaire intégrant une variable catégorielle qui indique si le patient est prioritaire.

Réponse :

Entrez votre texte ici

df1 <- df1 %>%
  mutate(prioritaire = ifelse(grepl("PRIO", Activity_DETAILS), "Oui", "Non"))

temps_systeme2 <- temps_systeme %>%
  dplyr::left_join(
    df1 %>%
      group_by(ID) %>%
      summarise(prioritaire = ifelse(any(grepl("PRIO", Activity_DETAILS)), "Oui", "Non")),
    by = "ID"
  )

df1 %>% filter(grepl("PRIO", Activity_DETAILS))

table(temps_systeme2$prioritaire)

Non Oui 
 45  19 

model2 <- lm(W ~ bloc_2h + prioritaire, data = temps_systeme2)
summary(model2)

Call:
lm(formula = W ~ bloc_2h + prioritaire, data = temps_systeme2)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5567 -0.8390 -0.1515  0.7885  3.7790 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)       0.93268    0.46253   2.016   0.0484 *
bloc_2h[10h;12h]  0.62588    0.52714   1.187   0.2399  
bloc_2h[12h;14h]  0.75690    0.58286   1.299   0.1992  
bloc_2h[14h;16h]  1.13320    0.56530   2.005   0.0497 *
bloc_2h[16h;18h]  1.00986    0.53958   1.872   0.0663 .
prioritaireOui    0.02496    0.31047   0.080   0.9362  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.126 on 58 degrees of freedom
Multiple R-squared:  0.08126,   Adjusted R-squared:  0.002062 
F-statistic: 1.026 on 5 and 58 DF,  p-value: 0.4109

# prédictions
temps_systeme2$pred_W2 <- predict(model2)

temps_systeme2 <- temps_systeme2 %>%
  mutate(prioritaire = factor(prioritaire, levels = c("Oui", "Non")))

temps_systeme2

ggplot(temps_systeme2, aes(x = bloc_2h, y = W, color = prioritaire)) +
  geom_jitter(width = 0.2, height = 0, alpha = 0.5) +
  geom_point(aes(y = pred_W2), shape = 17, size = 3) +
  labs(
    title = "Durée de séjour selon bloc 2h et priorité",
    x = "Bloc horaire d'arrivée",
    y = "Durée de séjour W (heures)",
    color = "Prioritaire"
  ) +
  theme_minimal()

NA
NA

Question 5) Modélisation de la durée de séjour par un taux de départ (/3)

Une manière alternative de “générer” un temps lié à un évènement est d’utiliser le taux de défaillance de sa distribution défini par : \[\mu(t) = \lim_{h \to +\infty} \frac{\mathbb{P}(X<t+h|X>t)}{h} = \lim_{h \to +\infty} \frac{\mathbb{P}(X<t+h) - \mathbb{P}(X<t)}{\mathbb{P}(X>t)h} = \frac{f(t)}{1-F(t)} = \frac{-\frac{dR(t)}{dt}}{R(t)} = -\frac{(ln(R(t))}{dt}\].

Dans le cas de la loi exponentiel, ce taux est constant car \(\frac{f(t)}{1-F(t)} = \lambda e^{-\lambda t} / (e^{-\lambda t}) = \lambda\), ce qui est une autre manière de voir que la la loi est sans mémoire.

Dans cette question, il vous est ainsi demandé :

• Exprimer (sous forme d’équation \(\mu(t)=...\)), calculer et tracer ce taux pour la distribution gamma pour les (3) distributions modéliser à la question 3 :

Réponse :

Pour une variable aléatoire \(X\) suivant une loi Gamma de paramètres shape \(k\) et rate \(\lambda\), le taux de défaillance est défini par :

\[ \mu(t) = \frac{f(t)}{1-F(t)} = \frac{\text{densité à t}}{\text{survie à t}}, \quad t>0 \]

où :

• \(f(t) = \dfrac{\lambda^k}{\Gamma(k)} t^{k-1} e^{-\lambda t}\) est la densité,
• \(F(t) = \int_0^t f(s) ds\) est la fonction de répartition,
• \(S(t) = 1 - F(t)\) est la fonction de survie.

Formules pour les distributions ajustées :

Patients prioritaires :

\[ \mu_{\text{prioritaires}}(t) = \dfrac{ \frac{1.857^{1.86}}{\Gamma(1.86)} t^{0.86} e^{-1.057 t} }{ 1 - \int_0^t \frac{1.857^{1.86}}{\Gamma(1.86)} s^{0.86} e^{-1.057 s} ds }, \quad t>0 \]

Patients non prioritaires :

\[ \mu_{\text{non-prioritaires}}(t) = \dfrac{ \frac{1.892^{3.257}}{\Gamma(3.257)} t^{2.257} e^{-1.892 t} }{ 1 - \int_0^t \frac{1.892^{3.257}}{\Gamma(3.257)} s^{2.257} e^{-1.892 s} ds }, \quad t>0 \]

Tous les patients :

\[ \mu_{\text{tous}}(t) = \dfrac{ \frac{1.5^{2.5}}{\Gamma(2.5)} t^{1.5} e^{-1.5 t} }{ 1 - \int_0^t \frac{1.5^{2.5}}{\Gamma(2.5)} s^{1.5} e^{-1.5 s} ds }, \quad t>0 \]

#_Entrez votre code R ici_
# Charger le package MASS pour fitdistr
library(MASS)

# Ajustement Gamma pour chaque groupe
fit_gamma_all       <- fitdistr(duree, "gamma")
fit_gamma_prio      <- fitdistr(duree_prio, "gamma")
fit_gamma_non_prio  <- fitdistr(duree_non_prio, "gamma")

# Récupération des paramètres shape et rate
shape_all  <- fit_gamma_all$estimate["shape"]
rate_all   <- fit_gamma_all$estimate["rate"]

shape_prio <- fit_gamma_prio$estimate["shape"]
rate_prio  <- fit_gamma_prio$estimate["rate"]

shape_non_prio <- fit_gamma_non_prio$estimate["shape"]
rate_non_prio  <- fit_gamma_non_prio$estimate["rate"]

# Définition de la fonction hazard pour la Gamma
hazard_gamma <- function(t, shape, rate) {
  f <- dgamma(t, shape=shape, rate=rate)   # densité
  S <- 1 - pgamma(t, shape=shape, rate=rate) # survie
  return(f / S)
}

# Vecteur de temps pour le calcul
t_vals <- seq(0.01, 10, by=0.01)  # éviter t=0 pour la Gamma

# Calcul des hazard pour chaque groupe
haz_all       <- hazard_gamma(t_vals, shape_all, rate_all)
haz_prio      <- hazard_gamma(t_vals, shape_prio, rate_prio)
haz_non_prio  <- hazard_gamma(t_vals, shape_non_prio, rate_non_prio)

# Tracé du taux de défaillance
plot(t_vals, haz_all, type="l", col="black", lwd=2,
     ylab="Taux de défaillance μ(t)", xlab="Temps t",
     main="Taux de défaillance - loi Gamma")
lines(t_vals, haz_prio, col="red", lwd=2)

lines(t_vals, haz_non_prio, col="blue", lwd=2)
legend("topright", legend=c("Tous", "Prioritaires", "Non prioritaires"),
       col=c("black","red","blue"), lwd=2)

NA
NA

• De construire et d’analyser un modèle (calcul de moyennes simples) qui estime ce taux de risque par tranche de 30min ([0;30min],]30min;60min]…) pour l’ensemble des patients, puis séparant selon le bloc d’arrivée de 2h [8h;10h], …, ]16h;18h], puis en séparant par patients prioritaires et non prioritaires

Entrez votre texte et formules ici

#_Entrez votre code R ici_

library(readxl)
library(dplyr)
library(lubridate)
library(MASS)

data_raw <- read_excel("Log_Patient_URO_12112015.xlsx")

# Conversion en format date-heure
data_raw <- data_raw %>%
  mutate(
    datetime_begin = ymd_hms(`Timestamp start`),
    datetime_end   = ymd_hms(`Timestamp end`)
  )

data_patient <- data_raw %>%
  group_by(ID) %>%
  summarise(
    arrival_time   = min(datetime_begin, na.rm = TRUE),
    departure_time = max(datetime_end, na.rm = TRUE),
    prioritaire = ifelse(any(grepl("PRIO", Activity_DETAILS)), "Oui", "Non"),
    .groups = "drop"
  ) %>%
  mutate(
    duree_sejour = as.numeric(
      difftime(departure_time, arrival_time, units = "hours")
    )
  )

# Nettoyage
data_patient <- data_patient %>%
  filter(!is.na(duree_sejour), duree_sejour > 0)

duree <- data_patient$duree_sejour

fit_norm  <- fitdistr(duree, "normal")
fit_exp   <- fitdistr(duree, "exponential")
fit_gamma <- fitdistr(duree, "gamma")
fit_weib  <- fitdistr(duree, "weibull")

# --- Calcul manuel de l’AIC
AIC_fitdistr <- function(fit) {
  k <- length(fit$estimate)
  -2 * fit$loglik + 2 * k
}

AIC_values <- data.frame(
  Distribution = c("Normale", "Exponentielle", "Gamma", "Weibull"),
  AIC = c(
    AIC_fitdistr(fit_norm),
    AIC_fitdistr(fit_exp),
    AIC_fitdistr(fit_gamma),
    AIC_fitdistr(fit_weib)
  )
)

AIC_values <- AIC_values[order(AIC_values$AIC), ]
print(AIC_values)

   Distribution      AIC
3         Gamma 186.0933
4       Weibull 190.3801
1       Normale 212.4402
2 Exponentielle 212.7157

breaks_30 <- seq(0, ceiling(max(duree)), by = 0.5)

hazard_par_tranche <- function(durees, breaks) {

  n <- length(breaks) - 1
  hazard <- numeric(n)

  for (i in 1:n) {

    debut <- breaks[i]
    fin   <- breaks[i + 1]

    a_risque <- sum(durees >= debut)
    sorties  <- sum(durees >= debut & durees < fin)

    hazard[i] <- ifelse(a_risque > 0, sorties / a_risque, NA)
  }

  data.frame(
    debut = breaks[-length(breaks)],
    fin   = breaks[-1],
    hazard = hazard
  )
}

hazard_global <- hazard_par_tranche(duree, breaks_30)

data_patient <- data_patient %>%
  mutate(
    heure_arrivee = hour(arrival_time),
    bloc_2h = cut(
      heure_arrivee,
      breaks = c(8, 10, 12, 14, 16, 18),
      right = FALSE,
      include.lowest = TRUE
    )
  )

hazard_par_bloc <- data_patient %>%
  group_by(bloc_2h) %>%
  group_map(~ hazard_par_tranche(.x$duree_sejour, breaks_30),
            .keep = TRUE)

hazard_priorite <- data_patient %>%
  group_by(prioritaire) %>%
  group_map(~ hazard_par_tranche(.x$duree_sejour, breaks_30),
            .keep = TRUE)

head(hazard_global)

hazard_par_bloc

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

hazard_priorite

[[1]]

[[2]]

library(ggplot2)

ggplot(hazard_global, aes(x = debut, y = hazard)) +
  geom_line(size = 1) +
  labs(
    title = "Taux de risque de sortie – Global",
    x = "Temps de séjour (heures)",
    y = "Taux de risque"
  ) +
  theme_minimal()

hazard_bloc_df <- bind_rows(
  hazard_par_bloc,
  .id = "bloc"
)

# Récupérer les vrais noms des blocs
hazard_bloc_df$bloc <- levels(data_patient$bloc_2h)[as.numeric(hazard_bloc_df$bloc)]

ggplot(hazard_bloc_df, aes(x = debut, y = hazard, color = bloc)) +
  geom_line(size = 1) +
  labs(
    title = "Taux de risque par bloc d’arrivée",
    x = "Temps de séjour (heures)",
    y = "Taux de risque",
    color = "Bloc d’arrivée"
  ) +
  theme_minimal()

hazard_prior_df <- bind_rows(
  hazard_priorite,
  .id = "prioritaire"
)

ggplot(hazard_prior_df, aes(x = debut, y = hazard, color = prioritaire)) +
  geom_line(size = 1) +
  labs(
    title = "Taux de risque – Patients prioritaires vs non prioritaires",
    x = "Temps de séjour (heures)",
    y = "Taux de risque",
    color = "Prioritaire"
  ) +
  theme_minimal()

Question 6) Calculs des niveaux d’occupation et validation des modèles (/5)

A partir des modèles précédant, il vous est demandé d’estimer, de visualiser et d’analyser le niveau d’occupation au cours de la journée en utilisant 30 échantillons de processus d’arrivée de patient (30 journées et pas 30 patients) Il vous est ensuite demandé de comparer la qualité des estimations des différents modèles à la réalité de manière visuelle et en utilisant une mesure d’erreur absolu moyenne et/ou quadratique (ou mis à la \(\sqrt()\)) et une mesure de biais sur l’occupation moyenne.

Réponse :

Entrez votre texte ici

#_Entrez votre code R ici_

# Fonction pour calculer l'occupation réelle
calculer_occupation_reelle <- function(data_patient, debut, fin, pas_temps = 0.1) {
  temps_seq <- seq(debut, fin, by = pas_temps * 3600)
  occupation <- sapply(temps_seq, function(t) {
    sum(data_patient$arrival_time <= t & data_patient$departure_time >= t)
  })
  data.frame(temps = temps_seq, occupation = occupation)
}

# Définir la période d'observation
date_jour <- as.Date("2015-11-12")
debut_journee <- ymd_hms(paste(date_jour, "08:00:00"))
fin_journee <- ymd_hms(paste(date_jour, "18:00:00"))

# Calculer l'occupation réelle
occupation_reelle <- calculer_occupation_reelle(data_patient, debut_journee, fin_journee)

# Taux d'arrivée par bloc de 2h
data_patient <- data_patient %>%
  mutate(
    heure_arrivee = hour(arrival_time),
    bloc_2h = cut(heure_arrivee, breaks = c(8, 10, 12, 14, 16, 18),
                  right = FALSE, include.lowest = TRUE)
  )

lambda_par_bloc <- data_patient %>%
  group_by(bloc_2h) %>%
  summarise(nb_arrivees = n(), lambda = nb_arrivees / 2)

# Paramètres des lois de durée
duree <- data_patient$duree_sejour
fit_exp <- fitdistr(duree, "exponential")
fit_gamma <- fitdistr(duree, "gamma")
fit_weib <- fitdistr(duree, "weibull")

# Fonction pour simuler une journée complète (arrivées + durées)
simuler_journee <- function(lambda_par_bloc, loi_duree, params_duree) {
  patients <- data.frame()
  
  for (i in 1:nrow(lambda_par_bloc)) {
    lambda <- lambda_par_bloc$lambda[i]
    debut_bloc <- c(8, 10, 12, 14, 16)[i]
    fin_bloc <- c(10, 12, 14, 16, 18)[i]
    
    # Nombre d'arrivées (Poisson)
    nb_arrivees <- rpois(1, lambda * 2)
    
    if (nb_arrivees > 0) {
      # Temps d'arrivée uniformes dans le bloc
      arrivees <- runif(nb_arrivees, min = debut_bloc, max = fin_bloc)
      
      # Durées selon la loi spécifiée
      if (loi_duree == "exp") {
        durees <- rexp(nb_arrivees, rate = 1/params_duree$mu)
      } else if (loi_duree == "gamma") {
        durees <- rgamma(nb_arrivees, shape = params_duree$shape, rate = params_duree$rate)
      } else if (loi_duree == "weibull") {
        durees <- rweibull(nb_arrivees, shape = params_duree$shape, scale = params_duree$scale)
      }
      
      patients <- rbind(patients, data.frame(
        arrivee = arrivees, duree = durees, depart = arrivees + durees
      ))
    }
  }
  return(patients)
}

# Fonction pour calculer l'occupation d'une simulation
calculer_occupation_simulation <- function(patients_sim, debut = 8, fin = 18, pas = 0.1) {
  temps_seq <- seq(debut, fin, by = pas)
  occupation <- sapply(temps_seq, function(t) {
    sum(patients_sim$arrivee <= t & patients_sim$depart >= t)
  })
  data.frame(temps = temps_seq, occupation = occupation)
}

# Lancer les 30 simulations pour chaque modèle
set.seed(123)
n_simulations <- 30

resultats_exp <- lapply(1:n_simulations, function(i) {
  calculer_occupation_simulation(
    simuler_journee(lambda_par_bloc, "exp", list(mu = 1/fit_exp$estimate["rate"]))
  )
})

resultats_gamma <- lapply(1:n_simulations, function(i) {
  calculer_occupation_simulation(
    simuler_journee(lambda_par_bloc, "gamma", 
                   list(shape = fit_gamma$estimate["shape"], rate = fit_gamma$estimate["rate"]))
  )
})

resultats_weib <- lapply(1:n_simulations, function(i) {
  calculer_occupation_simulation(
    simuler_journee(lambda_par_bloc, "weibull", 
                   list(shape = fit_weib$estimate["shape"], scale = fit_weib$estimate["scale"]))
  )
})

# Agrégation des résultats (moyenne et IC 90%)
agreger_simulations <- function(liste_resultats) {
  temps <- liste_resultats[[1]]$temps
  mat_occupation <- sapply(liste_resultats, function(x) x$occupation)
  data.frame(
    temps = temps,
    occupation_moy = rowMeans(mat_occupation),
    occupation_q05 = apply(mat_occupation, 1, quantile, probs = 0.05),
    occupation_q95 = apply(mat_occupation, 1, quantile, probs = 0.95)
  )
}

occupation_exp_agg <- agreger_simulations(resultats_exp)
occupation_gamma_agg <- agreger_simulations(resultats_gamma)
occupation_weib_agg <- agreger_simulations(resultats_weib)

occupation_reelle_plot <- occupation_reelle %>%
  mutate(temps_heure = as.numeric(difftime(temps, debut_journee, units = "hours")) + 8)

ggplot() +
  geom_line(data = occupation_reelle_plot, aes(x = temps_heure, y = occupation), 
            color = "black", linewidth = 1.2) +
  geom_line(data = occupation_exp_agg, aes(x = temps, y = occupation_moy), 
            color = "red", linewidth = 0.8, linetype = "dashed") +
  geom_ribbon(data = occupation_exp_agg, aes(x = temps, ymin = occupation_q05, ymax = occupation_q95),
              fill = "red", alpha = 0.1) +
  geom_line(data = occupation_gamma_agg, aes(x = temps, y = occupation_moy), 
            color = "blue", linewidth = 0.8, linetype = "dashed") +
  geom_ribbon(data = occupation_gamma_agg, aes(x = temps, ymin = occupation_q05, ymax = occupation_q95),
              fill = "blue", alpha = 0.1) +
  geom_line(data = occupation_weib_agg, aes(x = temps, y = occupation_moy), 
            color = "green", linewidth = 0.8, linetype = "dashed") +
  geom_ribbon(data = occupation_weib_agg, aes(x = temps, ymin = occupation_q05, ymax = occupation_q95),
              fill = "green", alpha = 0.1) +
  labs(title = "Occupation : Réel vs Modèles simulés (30 journées, IC 90%)",
       subtitle = "Noir = Réel | Rouge = Exp | Bleu = Gamma | Vert = Weibull",
       x = "Heure", y = "Nombre de patients") +
  scale_x_continuous(breaks = seq(8, 18, by = 2)) +
  theme_minimal()

calculer_metriques <- function(occupation_sim_agg, occupation_reelle_plot) {
  occupation_reelle_interp <- approx(
    x = occupation_reelle_plot$temps_heure,
    y = occupation_reelle_plot$occupation,
    xout = occupation_sim_agg$temps
  )$y
  
  mae <- mean(abs(occupation_sim_agg$occupation_moy - occupation_reelle_interp), na.rm = TRUE)
  rmse <- sqrt(mean((occupation_sim_agg$occupation_moy - occupation_reelle_interp)^2, na.rm = TRUE))
  biais <- mean(occupation_sim_agg$occupation_moy, na.rm = TRUE) - 
           mean(occupation_reelle_interp, na.rm = TRUE)
  
  c(MAE = mae, RMSE = rmse, Biais = biais)
}

tableau_metriques <- data.frame(
  Modèle = c("Exponentielle", "Gamma", "Weibull"),
  MAE = c(calculer_metriques(occupation_exp_agg, occupation_reelle_plot)["MAE"],
          calculer_metriques(occupation_gamma_agg, occupation_reelle_plot)["MAE"],
          calculer_metriques(occupation_weib_agg, occupation_reelle_plot)["MAE"]),
  RMSE = c(calculer_metriques(occupation_exp_agg, occupation_reelle_plot)["RMSE"],
           calculer_metriques(occupation_gamma_agg, occupation_reelle_plot)["RMSE"],
           calculer_metriques(occupation_weib_agg, occupation_reelle_plot)["RMSE"]),
  Biais = c(calculer_metriques(occupation_exp_agg, occupation_reelle_plot)["Biais"],
            calculer_metriques(occupation_gamma_agg, occupation_reelle_plot)["Biais"],
            calculer_metriques(occupation_weib_agg, occupation_reelle_plot)["Biais"])
)

print(tableau_metriques, row.names = FALSE, digits = 3)

        Modèle  MAE RMSE  Biais
 Exponentielle 1.54 1.92 -0.708
         Gamma 1.54 1.81  0.187
       Weibull 1.55 1.89  0.189

Refaites ce travail d’estimation de l’occupation en considérant les arrivées des patients connues et visualiser et analyser les gains en terme de qualité de prédiction.

Réponse :

Entrez votre texte ici

#_Entrez votre code R ici_

# Fonction de simulation avec arrivées connues
simuler_avec_arrivees_connues <- function(data_patient, loi_duree, params_duree, n_sim = 30) {
  # Extraire les heures d'arrivée réelles
  arrivees_reelles <- hour(data_patient$arrival_time) + minute(data_patient$arrival_time)/60
  n_patients <- nrow(data_patient)
  
  # Simuler n_sim journées
  lapply(1:n_sim, function(i) {
    # Générer les durées selon le modèle
    if (loi_duree == "exp") {
      durees <- rexp(n_patients, rate = 1/params_duree$mu)
    } else if (loi_duree == "gamma") {
      durees <- rgamma(n_patients, shape = params_duree$shape, rate = params_duree$rate)
    } else if (loi_duree == "weibull") {
      durees <- rweibull(n_patients, shape = params_duree$shape, scale = params_duree$scale)
    }
    
    # Calculer l'occupation
    calculer_occupation_simulation(
      data.frame(arrivee = arrivees_reelles, duree = durees, depart = arrivees_reelles + durees)
    )
  })
}

# Simulations pour chaque modèle
set.seed(456)
resultats_exp_connues <- simuler_avec_arrivees_connues(
  data_patient, "exp", list(mu = 1/fit_exp$estimate["rate"]), n_simulations
)

resultats_gamma_connues <- simuler_avec_arrivees_connues(
  data_patient, "gamma", 
  list(shape = fit_gamma$estimate["shape"], rate = fit_gamma$estimate["rate"]), 
  n_simulations
)

resultats_weib_connues <- simuler_avec_arrivees_connues(
  data_patient, "weibull", 
  list(shape = fit_weib$estimate["shape"], scale = fit_weib$estimate["scale"]), 
  n_simulations
)

# Agrégation
occupation_exp_connues_agg <- agreger_simulations(resultats_exp_connues)
occupation_gamma_connues_agg <- agreger_simulations(resultats_gamma_connues)
occupation_weib_connues_agg <- agreger_simulations(resultats_weib_connues)

# ==============================================================================
# VISUALISATION
# ==============================================================================

ggplot() +
  geom_line(data = occupation_reelle_plot, aes(x = temps_heure, y = occupation), 
            color = "black", linewidth = 1.2) +
  geom_line(data = occupation_exp_connues_agg, aes(x = temps, y = occupation_moy), 
            color = "red", linewidth = 0.8, linetype = "dashed") +
  geom_ribbon(data = occupation_exp_connues_agg, 
              aes(x = temps, ymin = occupation_q05, ymax = occupation_q95),
              fill = "red", alpha = 0.15) +
  geom_line(data = occupation_gamma_connues_agg, aes(x = temps, y = occupation_moy), 
            color = "blue", linewidth = 0.8, linetype = "dashed") +
  geom_ribbon(data = occupation_gamma_connues_agg, 
              aes(x = temps, ymin = occupation_q05, ymax = occupation_q95),
              fill = "blue", alpha = 0.15) +
  geom_line(data = occupation_weib_connues_agg, aes(x = temps, y = occupation_moy), 
            color = "green", linewidth = 0.8, linetype = "dashed") +
  geom_ribbon(data = occupation_weib_connues_agg, 
              aes(x = temps, ymin = occupation_q05, ymax = occupation_q95),
              fill = "green", alpha = 0.15) +
  labs(title = "Occupation avec ARRIVÉES CONNUES (30 simulations, IC 90%)",
       subtitle = "Noir = Réel | Rouge = Exp | Bleu = Gamma | Vert = Weibull",
       x = "Heure", y = "Nombre de patients") +
  scale_x_continuous(breaks = seq(8, 18, by = 2)) +
  theme_minimal()

# ==============================================================================
# MÉTRIQUES ET COMPARAISON
# ==============================================================================

# Métriques avec arrivées connues
metriques_exp_connues <- calculer_metriques(occupation_exp_connues_agg, occupation_reelle_plot)
metriques_gamma_connues <- calculer_metriques(occupation_gamma_connues_agg, occupation_reelle_plot)
metriques_weib_connues <- calculer_metriques(occupation_weib_connues_agg, occupation_reelle_plot)

# Tableau comparatif AVANT/APRÈS
comparaison <- data.frame(
  Modèle = rep(c("Exponentielle", "Gamma", "Weibull"), 2),
  Approche = c(rep("Arrivées simulées", 3), rep("Arrivées connues", 3)),
  MAE = c(
    calculer_metriques(occupation_exp_agg, occupation_reelle_plot)["MAE"],
    calculer_metriques(occupation_gamma_agg, occupation_reelle_plot)["MAE"],
    calculer_metriques(occupation_weib_agg, occupation_reelle_plot)["MAE"],
    metriques_exp_connues["MAE"],
    metriques_gamma_connues["MAE"],
    metriques_weib_connues["MAE"]
  ),
  RMSE = c(
    calculer_metriques(occupation_exp_agg, occupation_reelle_plot)["RMSE"],
    calculer_metriques(occupation_gamma_agg, occupation_reelle_plot)["RMSE"],
    calculer_metriques(occupation_weib_agg, occupation_reelle_plot)["RMSE"],
    metriques_exp_connues["RMSE"],
    metriques_gamma_connues["RMSE"],
    metriques_weib_connues["RMSE"]
  )
)

print(comparaison, row.names = FALSE, digits = 3)

        Modèle          Approche  MAE RMSE
 Exponentielle Arrivées simulées 1.54 1.92
         Gamma Arrivées simulées 1.54 1.81
       Weibull Arrivées simulées 1.55 1.89
 Exponentielle  Arrivées connues 1.03 1.30
         Gamma  Arrivées connues 1.39 1.81
       Weibull  Arrivées connues 1.39 1.78

# Calcul des gains
gains <- data.frame(
  Modèle = c("Exponentielle", "Gamma", "Weibull"),
  Reduction_MAE_pct = c(
    (calculer_metriques(occupation_exp_agg, occupation_reelle_plot)["MAE"] - 
     metriques_exp_connues["MAE"]) / 
     calculer_metriques(occupation_exp_agg, occupation_reelle_plot)["MAE"] * 100,
    
    (calculer_metriques(occupation_gamma_agg, occupation_reelle_plot)["MAE"] - 
     metriques_gamma_connues["MAE"]) / 
     calculer_metriques(occupation_gamma_agg, occupation_reelle_plot)["MAE"] * 100,
    
    (calculer_metriques(occupation_weib_agg, occupation_reelle_plot)["MAE"] - 
     metriques_weib_connues["MAE"]) / 
     calculer_metriques(occupation_weib_agg, occupation_reelle_plot)["MAE"] * 100
  ),
  Reduction_RMSE_pct = c(
    (calculer_metriques(occupation_exp_agg, occupation_reelle_plot)["RMSE"] - 
     metriques_exp_connues["RMSE"]) / 
     calculer_metriques(occupation_exp_agg, occupation_reelle_plot)["RMSE"] * 100,
    
    (calculer_metriques(occupation_gamma_agg, occupation_reelle_plot)["RMSE"] - 
     metriques_gamma_connues["RMSE"]) / 
     calculer_metriques(occupation_gamma_agg, occupation_reelle_plot)["RMSE"] * 100,
    
    (calculer_metriques(occupation_weib_agg, occupation_reelle_plot)["RMSE"] - 
     metriques_weib_connues["RMSE"]) / 
     calculer_metriques(occupation_weib_agg, occupation_reelle_plot)["RMSE"] * 100
  )
)

cat("\n=== GAINS EN QUALITÉ DE PRÉDICTION ===\n")

=== GAINS EN QUALITÉ DE PRÉDICTION ===

print(gains, row.names = FALSE, digits = 1)

        Modèle Reduction_MAE_pct Reduction_RMSE_pct
 Exponentielle                33              32.04
         Gamma                10              -0.08
       Weibull                10               5.91

Notes : Ici, on ne sépare pas les données en données d'entraînement et de test pour des raisons pratiques de quantité de données mais c'est un point à tenir en compte dans une vrai validation de modèle.

Question bonus (/1) :

Quelles sont les notations de Kendall des différents modèles de “file d’attente” que vous avez implémentés ? (Justifiez)

Note :  On est ici sur un cas spécial le modèle de file d'attente est sans file d'attente